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We develop finite temperature theory for a trapped dipolar Bose gas including thermal exchange interactions. 
Previous treatments neglected these, difficult to compute, terms. We present a methodology for numerically 
evaluating the thermal exchange contributions, making use of cylindrical symmetry. We then investigate prop- 
erties of the dipolar gas, including calculating the excitation spectrum over the full range of trap anisotropy. We 
evaluate the contributions due to thermal exchange noting that, under some regimes, these effects can be at least 
as significant as the direct interaction. We therefore provide guidance as to when these cumbersome terms can 
be neglected and when care should be exercised regarding their omission. 
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I. INTRODUCTION 

The creation of a Bose-Einstein condensate (BEC) of 52 Cr 
atoms in 2005 [ 1 1 has sparked many experimental and theo- 
retical investigations of the effect of long range dipole inter- 
actions on BECs. More recently, both 164 Dy flf) and 168 Er 
have been Bose-condensed. These three atomic species have 
large magnetic dipole moments compared to alkali atoms. 
There has also been progress towards a condensed dipolar 
Bose gas of molecules [4 J which gives the potential for much 
larger dipole moments. A recent experiment J3) has ex- 
plored the effect of dipole interactions on collective excita- 
tions. Theoretical investigation of the effect of dipolar interac- 
tions has often been focused on zero temperature fl6j-[9] . These 
works have studied excitations of the dipolar condensate as 
well as its ground state properties. Previous approaches to 
studying finite temperature effects include path integral Monte 
Carlo simulations IfTOHTTI and Hartree-Fock-Bogoliubov the- 
ory nn. 

We usually consider atoms in a BEC to interact via a van der 
Waals potential which drops off like 1/r 6 at large distances. 
It can be shown that for such a potential, the low energy scat- 
tering is dominated by the s-wave scattering [ 13 1. At temper- 
atures low enough for a BEC to form, only low energy scat- 
tering occurs so we can safely replace the actual interaction 
potential with a pseudopotential which reproduces the s-wave 
scattering. This pseudopotential is usually a delta function 
with the same s-wave scattering length. We cannot always 
assume that interactions between particles in a BEC are short- 
range only however. Chromium, dysprosium and erbium have 
significant magnetic dipole moments {6/j.b for 52 Cr Qj, 10/is 
for 164 Dy and l\i B for I68 Er 0, where fi B is the Bohr 
magneton). The dipolar interaction potential drops off like 
1/r 3 for large distances. In this case, all of the higher order 
partial waves contribute equally to the scattering at low en- 
ergy [13]. We therefore cannot replace the true potential with 
a short-range pseudopotential. When atoms with significant 
magnetic dipole moments are Bose-condensed we must take 
into account the long-range nature of the dipolar interaction 
between atoms and this gives rise to new and interesting be- 
haviour. 



Yi and You fl4l . and Goral et al. fl31 were the first to cal- 
culate the zero temperature density profiles for dipolar BECs 
using a dipolar Gross-Pitaevskii equation (GPE). Yi and You 
looked at the case of a pancake trap with A = y/8. They found 
that increasing the strength of dipolar interaction caused the 
condensate to contract in the radial direction and expand in the 
axial direction while the peak density increased. Meanwhile, 
Goral et al. compared the density profile of the dipolar BEC to 
an equivalent BEC with only contact interactions. They found 
that for a pancake trap, the dipolar BEC was larger in the ra- 
dial direction. For a cigar trap, the condensate was smaller in 
both the radial and axial direction. These studies showed that 
the anisotropic nature of the dipole interaction had significant 
effects on the shape of a condensate. The dipolar GPE was 
also used in the work of Santos et al. to study the stability 
of a dipolar condensate with no contact interactions. 

The natural next direction is to look at the excitation spec- 
trum of the system. Again, this problem was first approached 
by Yi and You fl6l . and Goral and Santos ifTTll . Both groups 
used a time-dependent, Gaussian variational ansatz in the GPE 
to calculate the energies of the three lowest excitations. Both 
groups also compare these results to numerical solutions of 
the time-dependent GPE. They do this by starting with the 
ground state wavefunction and applying a time-varying poten- 
tial to excite oscillations. Fourier analysis can then be applied 
to the changing condensate width to determine the oscillation 
frequency components. Both groups find good agreement be- 
tween variational and numerical results when the dipolar in- 
teraction is relatively weak, however, Goral and Santos find 
that this agreement breaks down when the dipolar interactions 
start to drive the condensate towards collapse. 

The standard approach to determining the excitation ener- 
gies of a BEC is to use Bogoliubov theory lITSTl . This approach 
was used by Ronen et al. to calculate the excitations of the 
dipolar BEC. The Bogoliubov method was not used by previ- 
ous studies lfT6l [T/l due to the numerical difficulty in solving 
the equations. Ronen et al. overcame this difficulty by util- 
ising both the cylindrical symmetry of the problem and the 
convenient form of the dipolar interaction term in momentum 
space. The authors calculated the excitation energies for vary- 
ing dipole strengths in a number of trap geometries. They also 
calculated the quantum depletion of the condensate. 



2 



In another article, Ronen et al. lfT9l examined the stability 
of the dipolar gas. They found that, contrary to a previous 
study |5), the condensate is unstable for any trap geometry if 
enough particles are added. They also found that for highly 
pancake traps with particular aspect ratios (e.g. A s» 7), the 
condensate develops a biconcave, or red blood cell-like, shape 
as it approaches instability. Ronen and Bohn |[T2l extended 
their method to condensates at non-zero temperature using 
Hartree-Fock-Bogoliubov theory with the Popov approxima- 
tion (HFBP) ll20ll2TI . The excitation energies and condensate 
fraction are calculated as a function of temperature. They also 
found that biconcave condensates still existed at finite temper- 
ature. 

In their finite temperature treatment, Ronen and Bohn ig- 
nore the effect of dipolar exchange from the thermal cloud on 
the condensate. The associated term in the HFBP equations 
is difficult to evaluate and was believed to be small. Here we 
include this term in order to determine whether it has any sig- 
nificant effects. In section [H] we will review the formalism 
of dipolar HFBP and the algorithm developed in earlier stud- 
ies [9, 12 1 to solve the associated equations for a cylindrically 
symmetric system. We will show how to calculate the thermal 
dipolar exchange term within this method. In sections [ill] and 
|PV| we calculate the density profiles and excitation energies 
for the system without including the thermal dipolar exchange 
term. In section [V] we then show the effect that the exchange 
term has on these properties. 



H. FORMALISM 

We describe the Bose condensed system using HFBP. This 
is the approach used by the authors of [12|. In that work, 
dipolar exchange interactions from the thermal cloud acting 
on the condensate were ignored due to the difficulty in eval- 
uating such terms. Here we calculate these terms in order to 
estimate how significant an effect they have on the properties 
of a condensate. The system is assumed to be in thermal equi- 
librium in the grand canonical ensemble. We work below the 
critical temperature and fix Nq atoms in the condensate. The 
chemical potential, /i, and the total number of atoms, N, are 
then determined by the temperature T. 

We examine a cylindrically symmetric system in a har- 
monic trapping potential, Vtr(x) = ^-(u> 2 p 2 + uj 2 z 2 ), where 
M is the atomic mass. We consider a gas of bosons that in- 
teract by both a long-range dipole-dipole interaction and a 
contact interaction characterized by V(r) = gS(r) + Vdd(r), 
where g — Anah 2 /M, with a the s-wave scattering length. 
We take the dipoles to be polarized along z giving a dipole 
interaction of 



Fdd(r) = 4^ lr|3 ' 



(1) 



where Cdd = MoMm f° r magnetic dipoles of strength \i m and 
d 2 /e for electric dipoles of strength d, and is the angle be- 
tween r and the z axis. 



The grand canonical Hamiltonian for this system is 

K = J dx*t( x ) (tf sp - M ) *(x) 

+ dxdx'&{x)&(x')V(x' -x)*(x')*(x) (2) 

where W(x) is the usual Bose field operator and _ff sp = 
(— fi 2 /2M)V 2 + Vtr(x) is the single particle Hamiltonian. 

The HFBP method consists of expanding the field opera- 
tor in a series of basis states and replacing the ground state 
operator with a c-number representing the condensate. We 
write this as ^(x) = v^o0o( x ) + ^( x )> where <fio(x) is the 
condensate wavefunction and i/>(x) is the fluctuation operator 
representing the other modes. When this is substituted into 
Eq. |2]) we can group the resulting terms based on the number 
of factors of ip(x) they contain. This gives terms up to quartic 
in ip(x), however, in the HFBP approach, the cubic and quar- 
tic terms are approximated using a mean-field factorisation 
giving terms which are quadratic or lower. We may therefore 
write the Hamiltonian as K = Kq + Ki + K^. By making 
the Hamiltonian stationary with respect to arbitrary first-order 
variation in ip(x), we obtain the generalised GPE for <fro(x) 



H sp + gn c (x) + 2gn(x) + $d(x) o (x) 

+ $ B [0 o (x)] =Wo( x ) (3) 



where g — Anh 2 a/M, n c (x) = -ZVo|0o( x )| 2 is the condensate 
density, and n(x) = (-0^( x )^( x )) is the non-condensate den- 
sity. The terms involving g come from the contact interaction, 
while $£i and $£; arise from the dipole-dipole interaction. 
The term $ o (x) represents the direct dipole interaction and 
has the form $d(x) = / dx.' Vd d ( x '- x )K( x ')+"-( x ')]- The 
other dipole term, <5>b[^)o( x )] represents the dipole exchange 
interaction between the condensate and non-condensate. It 
may be calculated from 

$ B [0o( x )]= y"dx'n(x',x)y dd (x'-x) ( /» (x') (4) 

where n(x',x) = (■0^( X ')V'( X )) is me non-condensate one 
body density matrix. 

In order to determine the excitations and non-condensate 
properties we turn to the second order contribution to the 
Hamiltonian, K-2- This may be diagonalised by applying a 
Bogoliubov transformation to the fluctuation operator, trans- 
forming to a new set of bosonic quasiparticle operators, 



Vi(x) 



E 

3 



v*(x.)a 



(5) 



and requiring Uj (x) and Vj (x) to satisfy the Bogoliubov -de 
Gennes (BdG) equations 



jCuj(x) — Mvj(x) = EjUj(x) 
C*vj(x) — M*Uj(x) = —EjVjix) 



(6) 
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where C = ho — p + M. ho is the operator on the left hand 
side of Eq. (j3j, i.e. ho4>o(x) = /i0o(x), while M is defined 
such that for an arbitrary function 4>(x), 

Mip(x) — gn c (x)ip(x)+ 

NoM*) I dx! o (x')Kid(x' - x)^(x') (7) 



Terms involving M represent exchange interaction with the 
condensate. The thermal one body density matrix can then be 
evaluated from 

= {K( x 'K(x) + Vj(x')v*(x)} n BE (E 3 ) 

3 

+v i (x')t$(x)} (8) 

where tjbe is the Bose distribution for the quasiparticles. The 
non-condensate density is simply h(x) = n(x, x). 

Following the method of l22l . to solve the BdG equations 
we first decouple equations |6} by introducing the new ampli- 
tudes ipf(x) = Uj(x) ± Vj(x). In terms of these amplitudes, 
the BdG equations are 



(h -(M + 2M)(h - m)^(x) = EU+(xl) 



(h - p) (ho - P + 2MU- (x) = EHt ( x ) 



3^3 
j^3 



(9a) 
(9b) 



Either of these equations can be used to determine the exci- 
tation energies, Ej, and original amplitudes Uj(x) and Vj(x). 
We choose to solve the first. We do this by expanding ^~(x) 
in the basis of eigenstates of the GPE ([3J, excluding the 
ground state, </>o(x), i.e. we write ^d"(x) = E a 4^a( x )' 

where (/io — p)cf) a (x) = e Q (/) Q ,(x). The eigenvalues, e Q are 
just the eigenvalues of the GPE relative to the chemical po- 
tential. This method has the advantage of ensuring that the 
excitations are orthogonal to the condensate. In this basis, 
Eq. (|9a|i is 



E 



(e a S 7a + 2M ia ) e a c* 



E?c? 



(10) 



the 2D Fourier transform in the x-y plane into a ID Hankel 
transform. The eigenstates of the GPE can be assumed to have 



the form /(x) 



p im<f> 



f(p, z) where to is an integer. Their 



Fourier transform, and any quantities with the same angular 
dependence, can be reduced to 

f(k p ,kt,k g ) = 2Tri- m e lmk * f dppf(p,z)J m (k pP ) (12) 



where J m (x) is the mth order Bessel function. We then define 
the mth order 2D Fourier-Hankel transform by 

Fm[f(p,z)] = 2tt J dppf(p,z)J m (k p p). (13) 

Previously, the thermal exchange term, in the GP and 
BdG equations has been ignored as it is computationally in- 
tensive to calculate and was not believed to be significant lfl2l . 
Here we calculate this term to determine the size of its effect. 
We can put Eq. Q in a similar form to ( fTT) 



®e [4>o (x)] 



1 



(2nf 



dke lkx l/ dd (k)/(x,k) (14) 



where /(x, k) = J dx' e~ jkx n(x' , x)0 o (x'). Unfortu- 
nately, Eq. ( fl4] l is not a Fourier transform due to the insep- 
arable x dependence in /(x, k), so the integral must be car- 
ried out for each grid position. To calculate excitations with 



to > we need to find $£[/(x)] for /(x) 



The Bogoliubov excitations correspond to distinct to values 



so we may write itj(x) 



Uj m (p,z) and likewise for 



Vj(x). Carrying out all of the angular integrals in Eq. ( |T4] > an- 
alytically, we get an expression for $ g which can be evaluated 
on a (p, z) grid 

* E [e im *f(p,z)] 

~im<f> 

~ ~2^T 



dk p k p j dk z e l " z Vdd(k p , k z )n m (p, z, k p , k z ), 

(15) 



where 



where M la = $ dx4>*(x)M(t> a (x). 

To solve Eq. Q we use the algorithm of |9|. With this 
method, the kinetic and dipolar interaction terms are calcu- 
lated using Fourier transforms. For example, the direct dipolar 
interaction is found from 



$i,(x) = J- 1 Vd d (k)Mk) + n(k)) 



(11) 



by using the convolution theorem. T~ l is the inverse Fourier 
transform and Vdd(k), n c (k) and n(k) are the Fourier trans- 
forms of the respective position space quantities. This method 
allows us to avoid the singular behaviour of Vdd (x) at the ori- 
gin. To improve accuracy, we truncated the interaction poten- 
tial in real space, using the method of ll23l where necessary. 

The cylindrical symmetry of the problem is used by per- 
forming the angular integrals analytically and thus converting 



(p, z, k pi k z ) ~ y ^ J-| m /_ m | [n m ' (p , z , p 7 z)]J\ m r_ m | (k p p) 



m'—0 

oc 



^ m' ~\-m 

[n* m ,(p', z', p, z)]J m , +m (k p p). 

m' — 1 

(16) 



and 



h m (p' , z' , p, z) = 

E i [ u im(p'. z')u jm (p, z) + v* m (p', z')v jm (p, z)] n B E(Ej) 



~ v j m (p', z')Vj m (p, z)} f(p', z'). 



(17) 



A number of schemes for choosing dimensionless parameters 
have been used in the literature. The energy and length units 
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FIG. 1 : Contour plots of condensate density for contact 
interactions only (solid) and contact and dipole interactions 
(dashed), with A = 1, T = hw p /k B , g = 1, N = 10 4 and 
D = 3 for the dipolar case. The density decreases 
monotonically from the origin, p, z are in units of aho- 



are usually chosen based on the harmonic trap frequency. Ex- 
cept where stated, a\ ia = y/K/MiJp. We also define 



1 



(2w 



(18) 



To characterise the strength of the dipolar and contact interac- 
tions we define the dimensionless parameters, 



D 



NMC dd 



and 



NMg _ AirNa 
ft 2 aho aho 



(19) 



(20) 



where N is the number of particles and a is the scattering 
length. 



III. DENSITY PROFILES 

As a basis for comparison, we first examine the solutions 
of the HFBP equations without including thermal dipolar ex- 
change. The first quantity we will look at is the condensate 
density. The presence of dipolar interactions alters the shape 
of the condensate as was seen in the earliest studies of dipo- 
lar condensates lfl4"l[T5l . The same behaviour occurs at finite 
temperature as we can see in Fig. [T] When there are only 
contact interactions present, a condensate in a spherical trap 
is also spherical. When there are dipolar interactions present, 
however, the condensate elongates in the axial direction. 

This behaviour can be explained by considering the 
anisotropy of the dipolar potential. Dipoles experience an at- 
tractive force when they are aligned head-to-tail rather than 
side-by-side. A density profile which is elongated along the 




(a) Cigar trap, A = 1/7 
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z 

(b) Pancake trap, A = 7 

FIG. 2: Condensate density with contact interactions only 
(solid) and with contact and dipolar interactions (dashed), (a) 
cigar-shaped trap with A = 1/7 and T = 0.5hui p /kB, (b) 
pancake trap with A = 7 and T = hu p /kB- All cases use 
N — 10 4 and g = 1. For the dipolar cases we used D = 3. 
The contour lines decrease in density from the centre 
outwards. Within each subfigure the solid and dashed contour 
lines are for the same density values, p, z are in units of aho- 



axial direction, so that more dipoles are aligned head-to-tail, 
will have a lower energy than one in which the dipoles experi- 
ence more repulsive forces. This is not a particularly intuitive 
result as one might expect that the repulsive interactions in 
the x-y plane would cause the condensate to expand radially. 
Since the dipoles are attractive along the axial direction how- 
ever, they may have a lower energy overall by aligning head- 
to-tail than they would by simply moving apart. The same 
behaviour was found at zero temperature by Yi and You lfl4ll 
for a slightly pancake trap. 

The effect of the dipolar interactions on the shape of the 
condensate is highly dependent on the aspect ratio of the trap 
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|[T5l . In Fig. [2] we have plotted the condensate density con- 
tours for a cigar trap with A = 1/7 and a pancake trap with 
A = 7. In the cigar-shaped trap, the condensate will be elon- 
gated in the axial direction. This means that the attractive 



part of the dipolar interaction will be dominant. In Fig. 2a 
we can see that the dipolar condensate drops off more steeply 
from the centre than the contact-only condensate. It also has 
a higher peak density. The attractive part of the dipolar inter- 
action is 'pulling' the condensate in on itself. In this case the 
dipolar interactions have little effect on the aspect ratio of the 
condensate. 



In the pancake-shaped trap in Fig. 2b b) we can see that the 



dipolar interaction affects both the drop-off in density and its 
aspect ratio. The condensate expands radially with the dipo- 
lar interactions and its density decreases more slowly from the 
centre than the contact-only case. This time the dipolar con- 
densate has a lower peak density than the contact-only con- 
densate. In the pancake trap, the repulsive part of the dipolar 
interaction is dominant and this explains these effects. Un- 
like in the spherically trapped situation, the condensate does 
expand radially due to the repulsion of side-by-side dipoles. 
In this case, the dipoles face a large energy cost for lining 
up head-to-tail due to the strong axial trapping. The lowest 
energy state is therefore when the dipoles spread out radially 
where there is only weak trapping. 

We can also look at the non-condensate density. This is 
given by setting x' = x in Eq. ([8) giving, 



n(x) 



{[k(x)| 2 + Mx)| 2 ]n BE (^) + k(x)| 2 }. (21) 



At zero temperature, tibe(E) = so the only contribu- 
tion to the non-condensate density is the quantum depletion, 
ri(x) = J2i |i!»(x)| 2 . This term generally requires the calcu- 
lation of very many modes to converge, however, it is usually 
very small [9|, and at finite temperature it is overwhelmed by 
the thermal population of the excited modes. The thermal part 
of the non-condensate converges more quickly as the Bose 
distribution factor decays exponentially with the mode energy 
once j3Ei 1. The fact that the quantum depletion may not 
be well converged is therefore not very important once there 
is an appreciable amount of thermal atoms. 

Two examples of non-condensate densities are shown in 
Fig. [3] The temperature used for these calculations is T — 
hujp/kB- This temperature is high enough that the contri- 
bution to the non-condensate density from thermal particles 
is much greater than that from the quantum depletion. It 
is also low enough that the energy cutoff of E cut — 6hui p 
does not significantly affect the result. Fig. [3a] shows the 
non-condensate density when there are only contact interac- 
tions present. The density is spherically symmetric, increas- 
ing from the centre of the trap to a maximum at p 2 + z 2 w 1 
and decreasing thereafter. This ring shaped thermal cloud 
also occurs when no interactions at all are included. This is 
simply due to the fact that the lowest modes above the con- 
densate mode have peaks away from the centre of the trap. 
This ring shaped cloud is enhanced by the presence of con- 
tact interactions as the thermal cloud experiences a potential 
of 2gn c (x) from the condensate which pushes it out. Fig. 3b 



a- 1 




1 1.5 

z 

(a) Contact interactions only 




(b) Contact and dipole interactions 

FIG. 3: Non-condensate density (a) without and (b) with 
dipolar interactions. Light is the highest density and dark is 

the lowest density. Parameters used are the same as for 
Fig.[T] Some of the contours are jagged due to the limited 
spatial resolution, p, z are in units of aho- 



shows the non-condensate density with dipolar interactions in- 
cluded. Compared with the non-dipolar case, the density is 
reduced in the z = plane of the trap and there are two sep- 
arate peaks at z/a^o m ±1. The two peaks arise due to the 
anisotropic dipolar interaction splitting up the ring from the 
contact interaction case. 

In this section we have shown the effects of dipolar interac- 
tions on the finite temperature BEC. The interactions can have 
a significant impact on the aspect ratio of the BEC, breaking 
the spherical symmetry when the trap is spherically symmet- 
ric. The effect of the interactions is strongly dependent on the 
aspect ratio of trap as this determines whether the attractive 
or repulsive part of the dipolar potential is dominant. In sec- 
tion [V] we will include the thermal dipolar exchange term and 
determine what further effect this has on the density of the 
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condensate. 



IV. EXCITATIONS 

We will now look at the excitation spectrum of the dipo- 
lar BEC. This is the set of eigenvalues, Ej, obtained from 
the BdG equations Previously, the zero temperature Bo- 
goliubov excitations have been studied by Ronen et al. Il9"l . 
Ronen and Bohn also looked at the finite temperature excita- 
tions without including thermal exchange lfl2l . They found 
that the excitation energies are not significantly affected by 
temperature except near the critical temperature or for highly 
non-spherical traps (e.g. A = 7). The temperatures we can 
accurately model are limited by the cutoff in energy we use 
for the Bogoliubov energies. Over this range of temperatures 
we also find very little change in excitation energies. We will 
therefore examine the zero-temperature spectrum as it is sig- 
nificantly easier to calculate and displays all of the same fea- 
tures as the small but finite temperature spectrum. We will 
also ignore the quantum depletion as its effect is very small 
and requires many modes to calculate accurately. 

First we examine how the excitations are affected by the 
aspect ratio of the trap. In order to concisely show the full 
range of trap aspect ratios, we follow the conventions used 
by Hutchinson and Zaremba 11241 . The aspect ratio is charac- 
terised by the parameter, 



Cutoff energy /fkjj p 



4 6 8 
0.82 0.92 0.99 



2u? 



(22) 



and energy and length scales are defined by the arithmetic 
mean as defined in Eq. ( [IS) , The aspect ratio, A, is related 

to fi by 



A 



'2/3 + 2 



(23) 



The parameter f3 can range from —1, corresponding to a 
one-dimensional cigar trap, to 2, corresponding to a two- 
dimensional pancake trap. The trap is spherical when (3 = 0. 
Because the trap potential is symmetric on reflection in the 
z = plane, the excitations can be split into those with even 
or odd amplitudes. The excitation energies are plotted against 
/3inFig.|4] 

The lowest m = 1 even excitation and lowest m = 
odd excitation represent the centre of mass or Kohn modes in 
the radial and axial directions respectively. The Kohn modes 
should oscillate at the trap frequencies, uj p for the radial 
mode and lj z for the axial mode, regardless of the interactions 
present. Their corresponding energies, in units of the mean 
trap energy, should therefore be given by uj p /uj a = \J\ — (3/2 
and uj z /u> a = ^/l + (3 for the radial and axial modes respec- 
tively. These analytical expressions have been plotted along 
with the calculated excitation energies in Fig. [4] There is ex- 
cellent agreement between the analytical and numerical re- 
sults indicating no problems with the model or numerics. 

The excitation energies for a non-interacting BEC are also 
plotted in Fig.|4]for comparison with the interacting case. The 



TABLE I: Ratio of contributions to the chemical potential 
from exchange and direct thermal dipolar interactions. Only 
modes which have an energy below the cutoff energy are 
included. The parameters used are N = 10 4 , X = V8, 
T = huj p /k B , D = 3.5 and g = 0. 



non-interacting excitations are simply calculated by setting 
both D and g to zero. The contact interaction used in the 
dipolar calculation is relatively small and the excitations pro- 
duced using only this interaction are only slightly different to 
the non-interacting case. The differences between the dipolar 
calculation and the non-interacting case are therefore primar- 
ily due to the dipole interaction. 

The dipolar interaction breaks the degeneracy between var- 
ious modes. This happens in two ways. Excitation energies 
corresponding to different values of m are split over a large 
range of the trap aspect ratio. Higher m excitations are pushed 
up relative to lower m excitations. For example, the lowest 
even m = 2 excitation is always degenerate with one of the 
even m = excitations in the non-interacting case. When 
dipolar interactions are included, the m = 2 excitation has a 
higher energy for most aspect ratios. As [3 approaches 2, and 
the trap becomes more pancake shaped, this energy splitting 
is reduced. This pattern breaks down for highly pancake traps 
both when the dipolar interaction is very strong [ 19] and at 
high temperatures fl2l . 

The dipolar interactions also break the degeneracy of modes 
of the same angular momentum at the aspect ratios where they 
would have crossed each other in the non-interacting case. 
This leads to avoided crossings between excitations of the 
same angular momentum. The energy differences are quite 
small making it difficult to see these avoided crossings in 
Fig. [4] 

In this section we have seen the effects of the dipolar inter- 
actions on the excitation spectrum of the Bose gas at zero tem- 
perature. The interaction breaks almost all of the degeneracies 
present in the low energy excitations for the ideal gas case, 
leaving only the degeneracy between +m and — m modes. In 
the next section we will examine the effect of the thermal ex- 
change term on these excitations at finite temperature. 



V. THERMAL EXCHANGE 

The dipolar thermal exchange term, ( p3| ), is the most dif- 
ficult to evaluate and we therefore attempt to minimise the 
number of times it is evaluated. We do this by first solving 
equations Q and ( |9a] > self-consistently without this term. This 
will be a good approximation to the full solution as long as the 
effect from dipolar thermal exchange is small. We then solve 
the full set of equations including dipolar thermal exchange 
using this solution as an initial guess. Assuming that the full 
solution is close to the initial guess, convergence should be 
achieved with relatively few self-consistency iterations. 
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(c) Even excitations, D = 0, g = 
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FIG. 4: (Color online) Zero temperature excitation energies in units of hui a (D and g use Oho = \Jh/Mu) a ). Different to states 
are represented by blue circles (to = 0), green squares (m = 1), red triangles (to = 2), cyan crosses (m = 3), and pink dots 
(to = 4). The solid black lines show the analytical values for the Kohn modes. 



To determine the size of the effect of the thermal exchange 
term we can examine the quantity $b[0o(x)]. This gives the 
interaction of the non-condensate on the condensate via ex- 
change. We cannot define an effective potential for the ex- 
change as we can with the direct dipolar interaction. 

To obtain a convenient measure of the importance of the 
thermal exchange term, we can look at its contribution to the 
chemical potential l25l . The chemical potential can be ob- 
tained from the GPE by taking the inner product with <f> (x) 
giving 

M = / 05( x )^s P <Mx) + [gn c (x.) + 2gh(x) 

+ $d(x)]|0 o (x)| 2 + 4>o (x) $ £ [0o (x)]dx . (24) 



The contribution of thermal exchange is therefore given by 

HE=f 0S(x)$ B [0o(x)]dx. (25) 

To determine the importance of this term we can compare it 
with the contribution from the direct dipolar interaction from 
the non-condensate. This is given by 

H 6 = J dx|0 o (x)| 2 J dx'Vd d (x'-x)n(x'). (26) 

We do not compare with the full direct dipolar interaction as 
this will be dominated by condensate-condensate interaction 
at low temperatures. 

We have examined the effect of the energy cutoff on the 
relative sizes of the direct and exchange thermal effects to de- 
termine how important different modes are to each effect. In 
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FIG. 5: (Color online) Contribution to to the chemical 
potential (in units of huj p ) from the direct (solid line) and 
exchange (dashed line) thermal dipolar interactions for a few 

iterations after thermal exchange is turned on. The trap 
aspect ratios used are A = v8 (blue crosses), A = 1 (green 
circles) and A = l/v8 (red squares). The other parameters 
used are N = 100, T = huj p /k B , D = 3.5 and g = 0. 




FIG. 6: (Color online) Condensate density cross-section 
along the z-axis extrapolated to p = 0. Parameters used are 
the same as in Fig. [5] with trap aspect ratios A = y/8 (blue), 

A = 1 (green) and A = 1/ v8 (red). The density before 
thermal dipolar exchange is included is given by the solid 
lines, while the dashed lines show the density after 
convergence with thermal dipolar exchange. 



Tab.|I] we show the ratio, \[ie/ M£>l> between the thermal and 
direct dipolar exchange contributions to the chemical poten- 
tial. As the cutoff energy is increased from four times to eight 
times the temperature, the relative importance of the exchange 
effect increases. This indicates that the higher energy, less 
occupied modes contribute more to the exchange effect than 
the direct effect. It is not clear that this trend will hold for 
higher energy modes which will become important at higher 
temperatures, however the HFBP calculations rapidly become 
infeasible as the energy cutoff is increased further. 

In order to present results for a system where a significant 
fraction of the atoms are in the thermal cloud, we assume a 
very small condensate with only 100 atoms. This allows us 
to use a relatively low temperature, and hence energy cutoff, 
while keeping 1-10% of atoms in the thermal cloud. This is 
enough for the thermal cloud to have a significant effect, but 
not so much as to require many iterations for self-consistency. 
We have performed calculations for three different trap aspect 
ratios: cigar-shaped, pancake-shaped and spherical. 

In Fig.[5]we show the contribution to the chemical potential 
from the direct and exchange thermal dipolar terms. We have 
plotted these for several iterations of the algorithm with the 
thermal exchange term included and we can see that there is 
reasonable convergence with only a few iterations. The contri- 
bution to the chemical potential from the direct thermal dipo- 
lar term is positive for a A = \/8 pancake trap as the repulsive 
part of the dipolar potential dominates in this trap. The con- 
tribution from the thermal exchange term is negative however, 
and of roughly the same magnitude. In a spherical trap the 
two terms both have a negative contribution of very similar 
magnitude. The size of both effects is reduced slightly com- 
pared to the pancake trap case. Finally, for the case of a cigar 



trap with A = 1 /V$>, the direct thermal dipolar contribution is 
strongly negative as the attractive part of the dipolar potential 
dominates. The thermal dipolar exchange term is again nega- 
tive and smaller in magnitude than the previous cases. In the 
cigar trap, the size of the exchange term is much smaller than 
the direct effect. 

In Fig. [6] we plot density cross-sections of the condensate 
for the same three cases. In each case, the central density 
is increased slightly with the inclusion of thermal dipolar ex- 
change. The effect is strongest for the pancake shaped trap, 
with the difference becoming imperceptible on the scale of 
the figure for the cigar trap. This matches the relative sizes of 
the thermal exchange effect on the chemical potential for the 
three trap aspect ratios. The importance of the thermal dipolar 
exchange term appears to be quite dependent on the trap as- 
pect ratio with the strongest effects in a pancake shaped trap. 
The thermal dipolar exchange term decreases the chemical po- 
tential and increases the central density of the condensate in 
each of the three cases shown. This shows that it is primarily 
an attractive interaction as opposed to the direct dipolar inter- 
action which can be attractive or repulsive depending on the 
aspect ratio of the trap. 

To show the effect that the thermal exchange interaction 
has on the excitation energies, we have calculated the shift in 
excitation energies when it is turned on. Tab. [Il] shows the 
excitation energies and their shifts for a spherical trap. The 
shifts in energy are very small and this is to be expected as 
there is little shift in the excitation energies at small but finite 
temperatures even when thermal dipolar exchange is excluded 
ifTZll . The two Kohn modes (with Ej s» 1) differ from the 
expected value slightly due to limitations of the HFBP method 
where the thermal cloud is static G6l . The results for the cigar 



9 



m Parity Ej AEj 

Odd 1.00 4.9 x 10~ 3 

Even 1.72 3.1 x 10" 4 

Even 2.15 1.2 x 10" 2 

1 Even 1.01 4.0 x 10" 3 

1 Odd 2.19 9.6 x 10" 3 

2 Even 2.39 8.7 x 10~ 3 



TABLE II: Shifts in the excitation energies due to thermal 
dipolar exchange for a spherical trap with the same 
parameters as in Fig. [5] The six lowest energy excitations are 
shown. The energies, Ej, are calculated without the thermal 
dipolar exchange term and AEj is the shift in these energies 
when this term is included. Energies are in units of huj p . 

and pancake shaped traps are similar with the effects again 
being largest for the pancake trap. 

The overriding conclusion however is that thermal ex- 
change can be at least as important as the direct interaction. 
This indicates care should always be exercised when neglect- 
ing thermal exchange. 

VI. CONCLUSIONS 

We have examined the effect of dipolar interactions on a 
harmonically trapped BEC. To begin with, we ignore the ther- 
mal dipolar exchange term. We have analysed the density pro- 
files of the condensate and thermal cloud. In a spherically 
symmetric trap we can see how the anisotropic dipolar inter- 
actions break the spherical symmetry of these density profiles. 
We calculate the zero-temperature Bogoliubov excitations of 
the system for a wide range of trap aspect ratios. Due to the 
limitation on the number of modes we can feasibly include 
in the calculations, we are limited to very low temperatures. 



We find that the excitation energies shift very little for these 
temperatures. The dipolar interactions break the degeneracy 
between a number of modes. 

We then include the thermal dipolar exchange interaction 
to determine its effect. This term cannot be expressed as an 
effective potential and it is much more computationally de- 
manding to calculate than the direct dipolar term. We have 
calculated the effect of the thermal dipolar exchange term on 
the chemical potential for different trap aspect ratios. We find 
that the term lowers the chemical potential and is most signif- 
icant for a pancake shaped trap. We also examine the effect 
on the density profile of the condensate. Again, the effect is 
largest for the pancake shaped trap and it tends to increase the 
central density. Finally, we determine what effect the thermal 
dipolar exchange term has on the excitation energies of the 
system. We find that the shift in the excitation energies due to 
thermal exchange is small. 

Overall, we find that the effect of the thermal dipolar ex- 
change can be as large as the thermal direct term at very low 
temperatures. These effects only become significant for very 
small condensates. Calculations for significantly higher tem- 
peratures become much more difficult as more modes need to 
be included. For the number of modes we have used (approx- 
imately fifty), calculations including exchange can take a few 
hundred hours. 
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